ylim.par = c(-1.0,1.0); xlim.par = c( 0.0,1.0)
no.excl_bene = 4; no.excl_cost = 2

rslt.file = paste(c(rslt.path,"r-conditional-identification-estimands.txt"),collapse = "")
mcb.local = array(dim = c(3,99,3))

for(i in 1:3){
	lower = (i - 1)*99 + 1; upper = lower + 98
	mcb.local[i,,] = as.matrix(read.table(rslt.file, na.strings = "."))[lower:upper, 2:4]
}
mcb.local = mcb.local/4

rslt.file = paste(c(rslt.path,"r-conditional-identification-support.txt"),collapse = "")
sup.local = array(dim = c(3,((no.excl_bene + 1) + (no.excl_cost + 1) + 1),2))

for(i in 1:3){
	lower = (i - 1)*((no.excl_bene + 1) + (no.excl_cost + 1) + 1) + 1; upper = i*((no.excl_bene + 1) + (no.excl_cost + 1) + 1) 
	sup.local[i,,] = as.matrix(read.table(rslt.file, na.strings = "."))[lower:upper,]
}
################################################################################
################################################################################
#illustrative graphs local identification (marginal benefits)
################################################################################
################################################################################
trace.pos = array(dim = c(max(no.excl_bene, no.excl_cost)))

for(h in 1:3){

if(h == 1){
	pdf(file 		 = paste(c(out.path,"g-conditional-identification-benefit-multiple.pdf"), collapse =""), width = 7, height = 5)
	shifters_name   = c("Short-Run Wages","Short-Run Unemployment","Tuition","Distance")
	shifters_legend = c("topleft","topright")
	no.local 	    = no.excl_bene 

	#position of label for exclusions
	trace.pos[1] = -0.40; trace.pos[2] = -0.55
	trace.pos[3] = -0.70; trace.pos[4] = -0.85

}

if(h == 2){
	pdf(file        = paste(c(out.path,"g-conditional-identification-cost-multiple.pdf"), collapse =""), width = 7, height = 5)
	shifters_name   = c("Long-Run Wages", "Long-Run Unemployment")
	shifters_legend = c("topleft","topright")
	no.local        = no.excl_cost 
	
	#position of label for exclusions
	trace.pos[1] = -0.4; trace.pos[2] = -0.55
}

if(h == 3){
	pdf(file        = paste(c(out.path,"g-conditional-identification-surplus-multiple.pdf"), collapse =""), width = 7, height = 5)
}

axes.par = FALSE
ylim.par = c(min(floor(mcb.local[c(1:3),,h])), max(ceiling(mcb.local[c(1,3),,h])))

for(g in c(1,3)){
	if(g != 1) axes.par = TRUE
	#----------------------------------------------------------------------------
	#marginal effects
	#----------------------------------------------------------------------------
	if(h == 1) k = 1
	if(h == 2) k = no.excl_bene + 2
	if(h == 3) k = no.excl_cost + no.excl_bene + 2 + 1

	plot(u.seq_double, mcb.local[g,x.idx_double,h], main = "",  pch = par.pch[g], xlab = "", ylab= "", col = "black", lwd = 1,
	     yaxs = "i", ylim = ylim.par, xlim = xlim.par, xaxs = "i", bty = "l");
	     lines(u.seq_single, mcb.local[g,,h])

	approx = approx(u.seq_double, mcb.local[g,x.idx_double,h], xout = seq(sup.local[g,k,1], sup.local[g,k,2], length   = 1000), method = "linear")
	par(new = TRUE)
	plot(approx$x, approx$y, main = "", pch = par.pch[g],xlab = "", 
		 ylab= "", col = "black", lwd = 1, yaxs = "i",
   	 ylim = ylim.par, xlim = c(0.0,1.0), xaxs = "i", bty = "l", axes = FALSE)
	#----------------------------------------------------------------------------
	#local identification
	#----------------------------------------------------------------------------
if(h < 3){
	for(i in 1:no.local){
 	 	k = k + 1; par(new = TRUE)
		x.par = seq(sup.local[g,k,1], sup.local[g,k,2], by = 0.02); y.par = rep(trace.pos[i], length(x.par))
		plot(x.par, y.par, main = "", pch = par.pch[i],xlab = "", ylab= "", col = "black", lwd = 1, 
			  yaxs = "i", ylim = ylim.par, xlim = xlim.par, xaxs = "i", bty = "l", axes = FALSE)
		lines(x.par, y.par)
		par(new = TRUE)
	}
}
par(new = TRUE)
}
#-------------------------------------------------------------------------------
#label
#-------------------------------------------------------------------------------
if(h == 1){; for(g in 1:2){
	if(g == 1){ upper = 2; lower = 1;}; if(g == 2){ upper = 4; lower = 3;};
	legend(shifters_legend[g], legend = shifters_name[lower:upper],
   	bty = "n", y.intersp = 1.5,  lty = 1  , pch = par.pch[c(lower,upper)],
   	col = "black", cex = 1.0)
}

	title(xlab = expression(u[S]), ylab = "Marginal Effects of Treatment")
	title(xlab = "         , p")
}

if(h == 2){ 
	lower = 1; upper = 2
	legend(shifters_legend[1], legend = shifters_name[lower:upper],
   	bty = "n", y.intersp = 1.5,  lty = 1  , pch = par.pch[c(lower,upper)],
   	col = "black", cex = 1.0)

	title(xlab = expression(u[S]), ylab = "Marginal Effects of Treatment")
	title(xlab = "         , p")
}

if(h == 3){
	title(xlab = expression(u[S]), ylab = "Marginal Effects of Treatment")
	title(xlab = "         , p")
}
	dev.off()
}
################################################################################
#Graph for Benefits and Cost						       									 
################################################################################
rslt.file = paste(c(rslt.path,"r-conditional-identification-single.txt"),collapse = "")

sup.idx = vector("numeric", length = 3)
sup.idx[1] = 1
sup.idx[2] = no.excl_bene + 1 + 1
sup.idx[3] = no.excl_bene + 1 + no.excl_cost + 1 + 1
	#----------------------------------------------------------------------------
	#benefit and cost graph
	#----------------------------------------------------------------------------
pdf(file=paste(c(out.path, "g-conditional-identification-cost-benefit-single.pdf"), collapse =""), width = 7, height = 5)

for(g in 1:2){
	#estimands
	axes.par = TRUE
	if(g == 2){; par(new = TRUE); axes.par = FALSE; }

   plot(u.seq_double, mcb.local[2,x.idx_double,g], main = "", 
        pch = par.pch[g],xlab = "", ylab= "", col = "black", lwd = 1,
        yaxs = "i", ylim = ylim.par, xlim = xlim.par, xaxs = "i",
        bty = "l", axes = axes.par);
	lines(u.seq_single, mcb.local[2,,g])
	#local identification
	approx = approx(u.seq_double, mcb.local[2,x.idx_double,g], xout = seq(sup.local[2,sup.idx[g],1], 
						 sup.local[2,sup.idx[g],2], length   = 1000),method = "linear")
	par(new = TRUE)
	plot(approx$x, approx$y, main = "",
	    pch = par.pch[g],xlab = "", ylab= "", col = "black", lwd = 1, yaxs = "i",
	    ylim = ylim.par, xlim = c(0.0,1.0), xaxs = "i", bty = "l", axes = FALSE)
}
title(xlab=expression(u[S]), ylab = "Marginal Treatment Effect")
      legend("bottom", legend = estimand.list[1:2], bty = "n", y.intersp = 1.5,
      lty = c(1,1,1)  , pch = par.pch[c(1,2)], col = "black", cex = 1.0)
par(new = FALSE); dev.off()
	#----------------------------------------------------------------------------
	#surplus graph
	#----------------------------------------------------------------------------
pdf(file=paste(c(out.path, "g-conditional-identification-surplus-single.pdf"), collapse =""), width = 7, height = 5)
plot(u.seq_double, mcb.local[2,x.idx_double,3], main = "", 
     pch = par.pch[1],xlab = "", ylab= "", col = "black", lwd = 1,
     yaxs = "i", ylim = ylim.par, xlim = xlim.par, xaxs = "i",
     bty = "l");
lines(u.seq_single,mcb.local[2,,3])

par(new = TRUE)
   
title(xlab=expression(u[S]), ylab = "Marginal Treatment Effect")
legend("bottom", legend = estimand.list[3], bty = "n", y.intersp = 1.5,
       lty = c(1,1,1)  , pch = par.pch[c(4)], col = "black", cex = 1.0)

approx = approx(u.seq_double,  mcb.local[2,x.idx_double,3], xout = seq(sup.local[2,sup.idx[3],1], sup.local[2,sup.idx[3],2], length   = 1000),method = "linear")
par(new = TRUE)
plot(approx$x, approx$y, main = "",
    pch = par.pch[1],xlab = "", ylab= "", col = "black", lwd = 1, yaxs = "i",
    ylim = ylim.par, xlim = c(0.0,1.0), xaxs = "i", bty = "l", axes = FALSE)
dev.off()




